% m_eval_1945_50_noise.m

% *********************************************************************************
% WE ADD SOME NOISE IN POPULATION AND EMPLOYMENT DISTRIBUTION IN 1945 
% RANDOMWLY ASSIGN NOISE (5%) FOR 100 TIMES AND USE THEM AS AN INITIAL DISTRIBUTION 
% *********************************************************************************
clear all; 
close all;
format shortG;
% ==================== LOAD FUNDAMENTAL AND DATA  =========================
load parameters.mat;                     % basic parameters 
load calib_productivity_amenity.mat;     % inverted fundamentals every 5 years from 1955
load calib_setup.mat;                    % Info data every 5 years from 1945
load commuting_cost.mat;  
load calib_auxiliary.mat;         % Info used in this calibration 

Rdata45 = Rdata_all(:,1); Rdata50 = Rdata_all(:,2); 
Ldata45 = Ldata_all(:,1); Ldata50 = Ldata_all(:,2); 
N = size(Ldata_all,1); 

% *************************************************************************
% SIMULATION SETUP 
% Generate noises on population and employment in 1945: for each simulation 
% randomly choose 50 locations for +10% population/employment and 50 locati
% ons for -10% population/employment and others are the same as in data   
% *************************************************************************
nSim=120; nNoise=50; 
rng(42); 
Noisevec=zeros(N,nSim); 

% Each column's rows are randomly ordered using sort of random keys
[~, idx] = sort(rand(N, nSim), 1);  % n x nSim; each column is a permutation
% Build linear indices for -1 and +1 placements
neg_lin = idx(1:nNoise,   :) + N*(0:nSim-1);  % k x nSims linear indices
pos_lin = idx(nNoise+1:2*nNoise,:) + N*(0:nSim-1);
Noisevec(neg_lin) = -0.05;
Noisevec(pos_lin) = 0.05;

all(sum(Noisevec == -0.05, 1) == 50);   % true: each column has 50 x -10%
all(sum(Noisevec == 0.05, 1) == 50);   % true: each column has 50 x +10%

Rdata45_noise = repmat(Rdata45,1,nSim).*(ones(N,nSim)+Noisevec);
Ldata45_noise = repmat(Ldata45,1,nSim).*(ones(N,nSim)+Noisevec); 

% =================== PARAMETERS ==========================================
alpha=alphaest;   
beta=betaest; 

% *************************************************************************
% ******** START SIMULATION ***********************************************
% *************************************************************************
R50_pmat = zeros(N,nSim); 
L50_pmat = zeros(N,nSim); 

for ss=1:nSim

R45_n = Rdata45_noise(:,ss); 
L45_n = Ldata45_noise(:,ss); 
Pop45_n = sum(R45_n); 
Emp45_n = sum(L45_n); 
Pop45 = min(Pop45_n,Emp45_n); 
R45_n = (R45_n./Pop45_n).*Pop45;
L45_n = (L45_n./Emp45_n).*Pop45;


% =========== COMPUTE LOWER BOUNDS for Theta 1950 & SET Theta 1950 ======== 
Rthetalb = (Rdata50 - R45_n < 0);   Rthetalb = (1 - Rdata50./R45_n).*Rthetalb;  Rthetalb = max(Rthetalb); 
Lthetalb = (Ldata50 - L45_n < 0);   Lthetalb = (1 - Ldata50./L45_n).*Lthetalb;  Lthetalb = max(Lthetalb);
thetalb = ceil(100.*max([Rthetalb; Lthetalb]))./100;

theta50 = 0.90;
theta55 = thetavec(1); 

%  ================ INVERSION OF (X,Y) FOR 1950 =============== 
options0 = optimset('Display','none','Algorithm','trust-region-dogleg','MaxFunEvals',50000000, 'MaxIter',1000000,'TolFun',1e-6,'TolX',1e-6);
X0=ones(N,1);  
Y0=ones(N,1);  
Rsyst = @(x)fReval2(x,K50,Q50,R45_n,Rdata50,L45_n,Ldata50,rho50,sigma50,theta50,mu,gammaL,gammaH);
Lsyst = @(y)fLeval2(y,K50,Q50,R45_n,Rdata50,L45_n,Ldata50,rho50,sigma50,theta50,mu,gammaL,gammaH);  
[x_fsolve, xfval_fsolve]=fsolve(Rsyst, X0, options0);   
[y_fsolve, yfval_fsolve]=fsolve(Lsyst, Y0, options0); 
X50=abs(x_fsolve); X50=X50./geomean(X50);
Y50=abs(y_fsolve); Y50=Y50./geomean(Y50);

% ************************************************
% *********** FORWARD LOOKING MODEL **************
% ************************************************
% =========== DECOMPOSE (X,Y) INTO FUNDAMENTALS AND AGGLOMERATION =========
% compute fundamental productivity and amenities in 1950
b50 = X50.*((Rdata50./Area).^(-beta)).*(X55.^(-rho*(1-theta55))); 
a50 = Y50.*((Ldata50./Area).^(-alpha)).*(Y55.^(-rho*(1-theta55))); 

% ========== COMPUTE STRUCTURAL RESIDUALS (FORWARD LOOKING) ================
l_ave_b = sum(log(bb),2)./size(bb,2); % average amenities over time (1955-75) 
l_ave_a = sum(log(aa),2)./size(aa,2); % average productivity over time (1955-75)

l_bb50 = log(b50) - l_ave_b;         % (log) structural error in amenities in 1950 (Forward Looking)
l_aa50 = log(a50) - l_ave_a;         % (log) structural error in productivity in 1950 (Forward Looking)

% Note: Since we have included the time-fixed fundamental for 1950, this
% is actually fundamentals + 1950 year fixed effects. 
l_b50_ave = sum(log(b50))./N; 
l_a50_ave = sum(log(a50))./N; 
bb_err50 = exp(log(b50)-l_b50_ave); 
aa_err50 = exp(log(a50)-l_a50_ave); 

XX50 = exp(l_ave_b).*((Rdata50./Area).^beta).*(X55.^(rho*(1-theta55)));   
YY50 = exp(l_ave_a).*((Ldata50./Area).^alpha).*(Y55.^(rho*(1-theta55))); 

% ********************************************************************
% ******** SIMULATE FORWARD LOOKING MODEL ****************************
% ******** WITHOUT STRUCTURAL RESIDUALS ******************************
% ** NOTE: HERE WE USE PROBABILITIES THAT IS BASED ON XX50, YY50 *****
% ********************************************************************
XX_temp = K50.*((repmat(Q50',N,1).^(-mu)).*repmat(XX50',N,1)).^(rho50/sigma50); 
XX_temp = XX_temp./repmat(sum(XX_temp,2),1,N); 

YY_temp = K50.*((repmat(Q50,1,N).^(-gammaH/gammaL)).*(repmat(YY50,1,N).^(1/gammaL))).^(rho50/sigma50); 
YY_temp = YY_temp./repmat(sum(YY_temp,1),N,1);

R50_p=Rdata50; 
L50_p=Ldata50; 
maxit=1e+6; 
ii=0; 
update=0.05;
while ii < maxit
    Ldif_p = L50_p-(1-theta50).*L45_n;
    Rdif_p = R50_p-(1-theta50).*R45_n; 
    
    R50_n = (1-theta50).*R45_n+sum(XX_temp.*repmat(Ldif_p,1,N),1)'; 
    L50_n = (1-theta50).*L45_n+sum(YY_temp.*repmat(Rdif_p',N,1),2);
    R50_p_r = round(R50_p*1e+6)./1e+6; R50_n_r = round(R50_n*1e+6)./1e+6; 
    L50_p_r = round(L50_p*1e+6)./1e+6; L50_n_r = round(L50_n*1e+6)./1e+6; 

    if min(R50_p_r == R50_n_r)==1 && min(L50_p_r == L50_n_r)==1
        break; 
    else
        ii=ii+1; 
        R50_p = (1-update).*R50_p+update.*R50_n; 
        L50_p = (1-update).*L50_p+update.*L50_n;
        if rem(ii,200)==1
        disp(['Largest Error is:',num2str(ii), ' ',num2str(max(abs(R50_p_r-R50_n_r))), ' ',num2str(max(abs(L50_p_r-L50_n_r)))]);
        end
    end
end 

R50_pmat(:,ss) = R50_p; 
L50_pmat(:,ss) = L50_p; 

disp(['Simulation Number is:',num2str(ss)]);

end 

% ==================  SAVE RESULTS ========================================
zeroCols = all(R50_pmat == 0, 1);    % columns that are entirely zeros
R50_pmat = R50_pmat(:, ~zeroCols);         % keep only non-zero columns
L50_pmat = L50_pmat(:, ~zeroCols);
nKept = size(R50_pmat, 2);

R50_sort = sort(R50_pmat,2); 
L50_sort = sort(L50_pmat,2); 
idxLow  = max(1, ceil(0.05 * nKept));   % often ceil for lower bound
idxHigh = min(nKept, floor(0.95 * nKept));
R50_low95  = R50_sort(:, idxLow);
R50_high95 = R50_sort(:, idxHigh);
L50_low95  = L50_sort(:, idxLow);
L50_high95 = L50_sort(:, idxHigh);

result=array2table([ID, X50, Y50, XX50, YY50, Rdata45, Rdata50, Ldata45,Ldata50, ...
    R50_low95, R50_high95, L50_low95, L50_high95, cbddist, Area, nKept.*ones(N,1)],'VariableNames',{'geocode','X50','Y50','XX50','YY50',...
    'pop45','pop50','emp45','emp50','pop50_low95','pop50_high95','emp50_low95','emp50_high95','cbddist', 'area','simtimes'}); 
writetable(result,'../../output/quant/output_evaluations_1945_50_noise.csv');

% ************************************************************************
% ************************************************************************
% **************************  FUNCTIONS **********************************
% ************************************************************************
% ************************************************************************
function[Reval] = fReval2(X,K,Q,DRpre,DRpost,DLpre,DLpost,rho,sigma,theta,mu,gammaL,gammaH)   

N = size(DRpre,1); 
X=abs(X); 
X=X./geomean(X);
xtemp=K.*((repmat(Q',N,1).^(-mu)).*repmat(X',N,1)).^(rho/sigma);
xtemp=xtemp./repmat(sum(xtemp,2),1,N); 
Rdif=DRpost-(1-theta).*DRpre; 
Ldif=DLpost-(1-theta).*DLpre; 
Reval = Rdif - sum(xtemp.*repmat(Ldif,1,N),1)'; 
Reval = abs(Reval); 

end

function[Leval] = fLeval2(Y,K,Q,DRpre,DRpost,DLpre,DLpost,rho,sigma,theta,mu,gammaL,gammaH) 

N = size(DRpre,1); 
Y=abs(Y);
Y=Y./geomean(Y);
ytemp= K.*((repmat(Q,1,N).^(-gammaH/gammaL)).*(repmat(Y,1,N).^(1/gammaL))).^(rho/sigma);
ytemp=ytemp./repmat(sum(ytemp,1),N,1); 
Rdif=DRpost-(1-theta).*DRpre; 
Ldif=DLpost-(1-theta).*DLpre; 
Leval = Ldif - sum(ytemp.*repmat(Rdif',N,1),2); 
Leval = abs(Leval); 

end
